% Scripts conducts a 3D rigid transform and cropping of ROI after a manual
% registration conducted in MeVisLab and after saving the metadata to a
% text-file.
% Optionally the calculation can be conducted with Matlab-routine, but it 
% is impossible to do so with big volume due to memory restrictions.
%
% Currently there is still an offset error (when running a rigid transform),
% that has to be found manually and additionally added to the metadata text
% file:
%--------------------------------------------------------------------------
% Date: 2014-01-14
% Author: Goran Lovric
% License: GPL 3
%--------------------------------------------------------------------------
clear;clc;close all;

src_list   = { '/sls/X02DA/Data20/e13657/temp_rigid/PSI-M2E-exv05', ...
               '/sls/X02DA/Data20/e13657/temp_rigid/PSI-M2E-exv10', ...
               '/sls/X02DA/Data20/e13657/temp_rigid/PSI-M2E-exv25'};
target_list = {'/sls/X02DA/Data20/e13657/temp_rigid/PSI-M2E-exv05crop3', ...
               '/sls/X02DA/Data20/e13657/temp_rigid/PSI-M2E-exv10crop3', ...
               '/sls/X02DA/Data20/e13657/temp_rigid/PSI-M2E-exv25crop3'};

metadata   = '/sls/X02DA/Data20/e13657/temp_rigid/metadata_final.txt';

with_matlab_routine = false;
tic
%--------------------------------------------------------------------------
% 0.) Load metadata
%--------------------------------------------------------------------------
fid = fopen(metadata);
tline = fgets(fid);
k = 1;
ii = 1;
while ischar(tline)
    A{k} = tline;
    if length(tline)>2
        if strcmp(tline(1:2),'==')
        headings(ii) = k;
        ii = ii+1;
        end
    end
    tline = fgets(fid);
    k = k + 1;
end
fclose(fid);

ind_3d = 1;
ind_transla = 1;
for kk = 1:length(headings)
        if strfind(A{headings(kk)},'Transform-matrix')
            for ii=1:4
                M_transfo(ii,:,ind_3d) = str2num(A{headings(kk)+ii});
            end
            ind_3d = ind_3d + 1;
        elseif strfind(A{headings(kk)},'Translation')
            T_transla(ind_transla,:) = str2num(A{headings(kk)+1});
            ind_transla = ind_transla + 1;
        elseif strfind(A{headings(kk)},'ROI')
            ROI(1,:) = str2num(A{headings(kk)+1});
            ROI(2,:) = str2num(A{headings(kk)+2});
        elseif strfind(A{headings(kk)},'WorldCoordinates')
            WorldCoordinates = str2num(A{headings(kk)+1});
        elseif strfind(A{headings(kk)},'RotationAxis')
            for ii=1:length(src_list)-1
                rotAxis(ii,:) = str2num(A{headings(kk)+ii});
            end
        elseif strfind(A{headings(kk)},'RotationCenter')
            for ii=1:length(src_list)-1
                rotCenter(ii,:) = str2num(A{headings(kk)+ii});
            end
        elseif strfind(A{headings(kk)},'EmpiricalCorrection')
            for ii=1:length(src_list)-1
                empCorrection(ii,:) = str2num(A{headings(kk)+ii});
            end
        end
end
if ~exist('empCorrection')
    empCorrection = zeros(length(src_list),3);
end

%--------------------------------------------------------------------------
% 1.) Load image stack
%--------------------------------------------------------------------------
for ll = 1:length(src_list)
src = src_list{ll};
target_dir = target_list{ll};
    
if exist(target_dir,'dir') ~= 7
    mkdir(target_dir);
end

files = dir([src  '/*.tif']);

zspace = '0000';

img = single(imread([src '/' files(1).name],'tif'));
img_sz = size(img);

I_mat = zeros(img_sz(1),img_sz(2),length(files),'single');

disp('')
disp('...loading data...')
for ind = 1:length(files)
    prog = ((ind/length(files))*100);
    if mod(prog,10) < 0.25;
        disp([num2str(round(prog)) ' % is done']);
    end
    zer=zspace(1:end-length(num2str(ind)));
	file=[src '/' files(ind).name];
    img = single(imread(file,'tif'));
	I_mat(:,:,ind) = img;
end

%--------------------------------------------------------------------------
% 2.) Rigid transform + ROI-cut
%--------------------------------------------------------------------------
x_range = [round(ROI(1,2)):round(ROI(2,2))];    % >> Mevislab Y
y_range = [round(ROI(1,1)):round(ROI(2,1))];    % >> Mevislab X
z_range = [round(ROI(2,3)):round(ROI(1,3))];    % >> Mevislab -Z

if ll > 1
    xx = M_transfo(2,4,ll-1)./WorldCoordinates;
    yy = M_transfo(1,4,ll-1)./WorldCoordinates;
    zz = M_transfo(3,4,ll-1)./WorldCoordinates;

    vec_rotAxis = -[rotAxis(ll-1,2) rotAxis(ll-1,1) rotAxis(ll-1,3)]./WorldCoordinates;
    vec_rotCentr= [0 0 0];%-[rotCenter(ll-1,2) rotCenter(ll-1,1) rotCenter(ll-1,3)]./WorldCoordinates;
    
    M_trans = AxelRot(rad2deg(rotAxis(ll-1,4)),vec_rotAxis,vec_rotCentr);
    
    if with_matlab_routine == true  % %>>>> with matlab routine
        tform = affine3d(M_trans.');
        I_mat = permute(I_mat,[2 1 3]);

        Rcb = imref3d(size(I_mat));
        Rout = Rcb;
        Rout.XWorldLimits = Rout.XWorldLimits -round(T_transla(ll-1,2)./WorldCoordinates);
        Rout.YWorldLimits = Rout.YWorldLimits -round(T_transla(ll-1,1)./WorldCoordinates); % numbers come from translation in RegistrationManual divided by world scaling
        Rout.ZWorldLimits = Rout.ZWorldLimits -round(T_transla(ll-1,3)./WorldCoordinates);  % Z again REVERTED

        new_img = imwarp(I_mat,tform,'OutputView',Rout);
        new_img = permute(new_img,[2 1 3]);
    else
        [new_img new_M] = affine(I_mat, M_trans);

        x_range = x_range-new_M(1,4)-round(T_transla(ll-1,2)./WorldCoordinates)+empCorrection(ll-1,1);
        y_range = y_range-new_M(2,4)-round(T_transla(ll-1,1)./WorldCoordinates)+empCorrection(ll-1,2);
        z_range = z_range-new_M(3,4)-round(T_transla(ll-1,3)./WorldCoordinates)+empCorrection(ll-1,3);
    end
else
    new_img = I_mat;
    new_M = zeros(4);
end
        
new_img = new_img(x_range,y_range,z_range);

%--------------------------------------------------------------------------
% 3.) Save images
%--------------------------------------------------------------------------
zspace = '0000';

disp('')
disp('...saving data...')
for ind = 1:size(new_img,3)
    prog = ((ind/size(new_img,3))*100);
    if mod(prog,10) < 0.25;
        disp([num2str(round(prog)) ' % is done']);
    end
    zer=zspace(1:end-length(num2str(ind)));
	file=sprintf('%s/%s%d.tif',target_dir,zer,ind);
    img = uint8( new_img(:,:,ind) ); % save as 8bit
	imwrite(img,file,'tif');
end

end

toc

% Release memory
clear new_img
clear I_mat